NASA TM X- 71918 




NASA TECHNICAL 

MEMORANDUM 


NASA TM X- 71918 

COPY NO. 12. 


THE COMPUTATION OF STEADY NOZZLE FLOW 
BY A TIME-DEPENDENT METHOD 

By Michael C, Cline 
April 197A 


(NASA-TM-X-71918) THE COMPUTATION OF N74-17013\ 

STEADY NOZZLE FLOW BY A TIME-'DEPENDENT 
METHOD (NASA) — p HC $4.00 CSCL 20D 

^ Unclas 

G3/12 29638 ) 


Reprodtieed by a i 

national technical 
information service 

U S Depert»an< ComtnBres 
Springfield, VA. 22151 


This informal documentation medium is used to provide acce! erated or 
special release of technical information to selected users. The contents 
may not meet NASA formal editing and publication standards, may be re- 
vised, or may be incorporated in another publication. 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 
LANGLEY RESEARCH CENTER, HAMPTON, VIRGINIA 23665 



2. Government Acceetion No. 


3. Recipient't Cataiog No. 


1. Report No. 

NASA TM X-71918 


15. Supplementary Notes 

Special technical information release, not planned for formal NASA publication 


16. Abstract 

The Steady flow in two-dimensional and axlsymmetric nozzles was computed 
using a time-dependent method* In this method the interior mesh points were 
computed using the MacComack finite-difference scheme, while a characteristic 
scheme was used to calculate the boundary mesh points. No explicit artificial 
viscosity term was included. The fluid was assumed to be a perfect gas. This 
method was used to compute the flow in a 45* - 15* conical, converging- diverging 
nozzle, a 15* conical, converging nozzle, and a 10* conical, plug nozzle. Good 
agreement between the numerical solution and experimental data was found. In 
contrast to previous time-dependent methods, the computational times were less 
than one minute on a CDC 6600 computer. 


4. Title and Subtitle 

The Computation of Steady Nozzle Flow By A Time- 
Dependent Method 

7. Author(s) 

Michael C. Cline 

9. Performing Organitation Name and Address 

NASA Langley Research Center 
Hampton, VA 23665 


12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D.C., 20546 &/ U.S. Atomic Energy 
Commission, Los Alamos, NM 87544 



6. Performing Organization Code 


8. Performing Organization Report No. 


IQ. Work Unit No. 


lit • ■ Crc^lfrraiii 


11. Contract or Grant No. 


13. Type of Report and Period Covered 

Technical Memorandum 


14. Sponsoring Agency Code 


17. Key Words (Suggested by Author(s)) (STAR category underlined) 

Aerodynamics 

Subsonic Flow Nozzle Flow 

Transonic Flow Channel Flow 

Supersonic Flow 
Hypersonic Flow 

18. Distribution Statement 

Unclassified > 

- Unlimited 


19. Security Qassif. (of this report) 
Unclassified 

L- — 

20. Security Classif. (of this page) 

Unclassified 

21. No. of Pages 

22. Price* 

NTIS ^,60 


S The National Technical Information Service, Springfield, Virginia 22151 
ST IF/NASA Scientific and Technical Information Facility. P.O. Box 33, College Park. MD 20740 


/ ) 











THE COMPUTATION OF STEADY NOZZLE FLOW 


BY A TIME- DEPENDENT METHOD* 

Michael C. Cline** 

University of California 
Los Alamos Scientific Laboratory 
Los Alamos, New Mexico 87544 

Abstract 

The steady flow in two-dimensional and axis ymme trie nozzles was 
computed using a time-dependent method. In this method the interior mesh 
points were computed using the MacCormack finite-difference scheme, while 
a characteristic scheme was used to calculate the boundary mesh points. 

No explicit artificial viscosity term was included. The fluid was assumed 
to be a perfect gas. This method was used to compute the flow in a 45® - 15® 
conicai., converging-diverging nozzle, a lO'' conical, converging nozzle, and 
a 10® conical, plug nozzle. Good agreement between the numerical solution 
and experimental data was found. In contrast to previous time- dependent 
methods, the computational times were less than one minute on a CDC 6600 
computer. 

*lhe initial phase of this research was accomplished while the author 
held a National Research Council Postdoctoral Resident Research Associate- 
ship supported by the NASA Langley Research Center, Hampton, Virginia. Tue 
final phase was performed under the auspices of the United States Atomic 
Energy Commission. 

Index categories: Subsonic and Transonic Flow; Supersonic and Hyper- 

sonic Flow; Nozzle and Channel Flow. 
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Introduction 


The equations of motion governing steady, Inviscld flow are of a 

mixed type, that is, hyperbolic In the supersonic region and elliptic in 

the subsonic region* These mathematical difficulties may be removed by 

using the so called time-dependent method* In this technique the flow is 

assumed to be unsteady or time-dependent. The governing equations are,. 

therefore, hyperbolic in both the subsonic and supersonic regions. The 

steady state solution may be obtained as the asymptotic solution for large 

time* This technique has been used to compute converging-diverging nozzle 
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flows by Prozan (as reported by Saunders and Cuffel, et al. ), Mlgdal, 

3 ^5 $ 

et al. , Wehofer and Moger , Laval , and Serra . This technique has also 

been lised to compute converging nozzle flows by Wehofer and Moger and 

Bro\>m and Ozcan^. While the results of the above calctilations are for the 

most part good, the computational times are rather large. In addition, 

although the computer program of Ref. 6 included a centerbody and those 

of Refs. 4 and 7 included the exhaust jet, none of the above codes is 

able to calculate both, that is, plug nozzles. Therefore, the object of 

this research was to develop a production type computer program capable 

of solving converging, converging- diverging, and plug two-dimensional 

nozzle flows in computational times of one minute or less on a CDC 6600 

computer. 

Literature Review 

The following is a discussion of the methods used in Refs. 1 through 
7. The first paragraph deals with the computation of the interior mesh 
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points while the next three paragraphs are concerned with the boundary 
tncsh points. 

Prozan , Wehofer and Moger, and Laval used variations of the two-step 
Lax-Wendroff scheme to compute the interior mesh points. Mlgdal, et al,, 
and Brown and Ozcan employed the original one-step Lax-Wendroff scheme, 
but with the equations of motion In non-conservation form, Serra applied 
the original Lax-Wendroff scheme with the equations of motion in conserva- 
tion form. In order to stabilize their schemes, Laval and Serra used 
®*^tificial viscosity terms in their difference equations. Wehofer and 
Moger reset the stagnation conditions along each streamline, reset the 
mass flow at each axial location, and smoothed the subsonic portion of 
the flow after each time step. 

In order to compute the nozzle inlet mesh points, Prozan^ assumed the 
inlet flow to be uniform, Wehofer and Moger assumed only that the pressure 
was radially uniform at the inlet. Migdal, et al. , and Brown and Ozcan 
mapped the inlet to minus infinity after Moretti , thus allowing the 
static conditions to be set equal to the stagnation conditions. Laval 
used extrapolation of the Interior mesh points to determine the inlet 
mesh points, while Serra employed a characteristic scheme. 

Prozan,^ Wehofer and Moger, Laval, and Brown and Ozcan used an extra- 
polation technique to compute the wall mesh points. Migdal, et al., em- 
ployed a characteristic scheme after Moretti to compute the wall mesh 
points, while Serra applied a reflection technique. In order for the con- 
verging nozzle problem to be properly posed, an exhaust jet calculation must 
be Included. Wehofer and Moger used an extrapolation procedure to compute 
Che exhaust jet boundary mesh points, while Brown and Ozcan employed a 

D 

characteristic scheme after Moretti . 



AH of the above authors used extrapolation to compute the exit 
mesh points when the flow was supersonic, since any errors Incurred 
would be swept out of the mesh* Serra employed a characteristic scheme 
when the exit flow was subsonic. 

Choice of Method 

The lengthy computational times associated with timer-dependent 
calculations are usually due to either inefficient numerical schemes or 
poor treatment of the boundaries resulting in the requirement for ex- 
cessively fine computational meshes. In the following paragraphs, a 
technique for the much more efficient calculation of the interior and 
boundary mesh points will be discussed. 

Ttie computation of steady flows by a time— dependent method differs 
from ordinary initial-value problems in that the initial data and much of 
the transient solution have a negligible effect on the final or steady 
solution. Therefore, accuracy is important only for the asymptotic state, 
and special attention can be given to intermediate efficiency In order 
that computational times are made reasonable* For this reason, the in- 
terior mesh points can be computed using a very efficient finite-difference 
scheme, as opposed to those less efficient finite- difference or character- 
istic schemes that achieve high accuracy at every step. 

In the class of finite- difference schemes, the two-step methods such 
as the MacCormack^^ and the two-step Lax-Wendroff schemes^^ are more ef- 
ficient than the original Laxp-Wendroff scheme^^, especially if the govern- 

12 

ing equations are in conservation form. Morettl showed that using the 
equations of motion in conservation form deceased efficiency and ease of 
programming while only slightly increasing the accuracy of shock calculations. 


-4 



The use of an explicit artificial viscosity term also decreases efficiency 

12 

and was shown to be physically unjustified by Moretti « In addition, such 
increases in the numerical dissipation can often destroy the weak shock 
structure of transonic flows. Therefore, the MacCormack scheme with the 
equations of motion in non-conservation form is used to calculate the in- 
terior mesh points. No explicit artificial viscosity term was included, 
although an implicit dissipation is present as an effect of truncation 
terms, which is sufficient to assure numerical stability for the present 
flows. 

Tlie boundary mesh points , while making up only a small part of the 

g 

total mesh points, are the most important to be treated accurately , be- 

Q 

cause of the flow-field sensitivity to precise boundary geometry. Moretti 

9 

and Abbett showed that reflection, extrapolation, and one-sided difference 
techniques for computing solid wall boundaries g.ive poor results and 
should be avoided. Therefore, the wall and centerbody mesh points are com- 
puted using a characteristic scheme. Likewise, the exhaust jet boundary 
mesh points are also calculated employing a characteristic scheme. 

In the case of the nozzle inlet mesh points, the use of extrapolation 
techniques and the assumption of one- dimensional flow presume the form of 
the solution and in many cases are physically unjustified. On the other 
hand, a characteristic scheme could be used to calculate the inlet mesh 
points. While the stagnation pressure and temperature are assumed to re- 
main constant at the inlet in a characteristic scheme, which is not necessar 
lly the case for unsteady flow, this assumption would appear to be valid 

o 

for the time-dependent calculation of steady flows. However, Moretti 
recommends mapping \he inlet to minus infinity, thus allowing the static 
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conditions to be set equal to the stagnation conditions. In theory, this 
would appear to be the best approach. However, it should be kept in tnlnd 
that the infinite physical plane taust be replaced by a finite computational 
plane. Also, this technique requires additional mesh points upstream of the 
nozzle inlet. It is not presently resolved as to whether the characteristic 
scheme approach us«d by Serra or the mapping to minus infinity approach 

g 

suggested by Moretti and employed by Migdal, et al, , and Brown and Ozcan 
Is the best technique. To reduce the total number of mesh points to be 
computed, a charactertlstlc scheme is used to compute the inlet mesh points. 

Extrapolation is used to compute the exit mesh points when the flow is 
supersonic and a characteristic scheme is employed when the flow is subsonic. 

Equations of Motion 

The appropriate non- conservation form of equations for two-dimensional, 
Inviscid, isentropic, rotational flow are 

+ up + VP + pu + pv + epv/y “ 0 (1) 

txyxy "^ 

+ vUy + p^/p = 0 (2) 

2 

^t *'■ “Px "^y “ ^ ^^t *'* '^^y^ ' ^ 

where p is the density, u is the axial velocity, v is the radial velocity, 
p is the pressure, a is the local speed of sound, t is the time, x and y 
are the axial and radial coordinates, and the subscripts denote partial 
differentiation. The symbol e is 0 for planar flow and 1 for 
axisymmetric flow. 
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1116 physical (x, y) plane is mapped into a rectangular computational 


plane (C, n) by the following coordinate transformation; 


C “ X ; n 


y-y (x) 
c 


y^(x.t) - y^(x) * 


T - t 


(5) 


where y^(x, t) denotes the nozzle wall and exhaust jet boundary radius as 
a function of x and t and y^(x) denotes the nozzle centerbody radius as a 
function of x. In the (c, n* t) coordinate system Eqs. (1) through (4) 
become 

P + up + vp + pu + pau + pSv + epv3/(y +n/3) “ 0 (6) 

^ s n si rj n c 


u + uu^ + vu^ + p^/p + op /p = 0 

T ^ n 4 n 


(7) 


V + uv^ + w + 6p /p = 0 
T c n n 


up^ + vp - a(p +UP +vp)«0 
T c n T c n 


( 8 ) 

(9) 


where 


1 

; a - - e 


y ** y 

■^w ^ c 


3x 
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6 = - n0 


3y, 


w 


3t 


( 10 ) 


V *= au + 3v + 6 (11) 

The fluid is assumed to be thermally and calorically perfect; that 
Is, a constant ratio of specific heats. 

Numerical Method 


The computational plane is divided up into five sets of mesh points. 
These five sets are the interior, inlet, exit, wall and centerbody, and 
exhaust jet boundary mesh points. The treatment of each of these set^v is 
given below. 
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Interior Mesh Points 


The interior mesh points are computed using the MacConnack scheme. 

This scheme is a second-order, noncentered, two-step, finite-difference 
scheme. Backward differences are used on the first step while forward 
differences are employed on the second step. The governing equations are 
left in non-conservation form. No explicit artificial viscosity term is 
used. Centerline mesh points are computed enforcing symmetry of the flow. 

A complete description of the method is given in Ref. 10. 

Inlet Mesh Points 

The inlet mesh points are computed using a second-order, reference- 
plane characteristic scheme. In this scheme the partial derivatives with 
respect to n are computed in the initial-value and solution surfaces using 
noncentered differencing as in the MacCormack scheme. These approximations 
to the derivatives with respect to n are then treated as forcing terms and 
the resulting system of equations is solved in the n “ constant reference 
planes using a two-independent variable, characteristic scheme. The scheme 
consists of a first step using backward differences in the initial-value 
plane and a second and final step using an average of the backward n dif- 
^®^Gt^ces in the Initial— value plane and forward n differences in the solution 
plane. The boundary condition is the specification of the stagnation temper- 
ature and stagnation pressure. The use of a reference-plane characteristic 
scheme requires an additional boundary condition. This additional boundary 
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condition Is the specification of the inlet flow angle, since In 
general the Inlet flow angle can be approximately detomdned from the 
nozzle geometry. These three quantities, along with the one Mach-line 
compatibility equation, are sufficient to determine the four dependent 
variables, u, v, p, and p. 

Exit Mesh Points 

Subsonic Flow . For subsonic flow a reference-plane characteristic 
scheme similar to the inlet scheme is used. The exit pressure is 
specified. This pressure, one Mach-line compatibility equation, and two 
streamline compatibility equations are sufficient to determine the four 
dependent variables. 

Supersonic Flow . For supersonic flow the exit mesh points are com- 
puted using linear extrapolation. Since any errors incurred here will be 
swept out of the mesh, such a procedure should be sufficient. 

Wall and Centerbody Mesh Points 

The wall and centerbody mesh points are also computed using a 
reference- plane characteristic scheme. In this scheme the derivatives with 
respect to ^ are approximated and the resulting system of equations is 
solved in the 5 « constant reference planes. The wall and centerbody 

contours, and therefore, their slopes are specified. The tangency conaltlon 
given by 

V “ u tan 0 + 3y^/St , (12) 

where 0 is the local wall or centerbody angle, one Mach-line compatibility 
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equation, and two streamline compatibility equations are sufficient to 
determine the dependent variables. 

Exhaust Jet Boundary Mesh Points 

The exhaust jet boundary mesh points are computed by the wall routine 

such that the pressure boundary condition 

p = p (13) 

*^ambxent 

is satisfied. This is accomplished by first assuming the shape of the 
jet boundary and then using the wall routine to calculate the pressure. 

Next, the jet boundary location is slightly changed and a second pressure 
is computed. By use of an interpolation-extrapolation procedure, a new 
jet boundary location is determined. This interpolation-extrapolation 
procedure is then repeated at each point until the jet boundary pressure 
and the ambient pressure agree to within some specified tolerance. 

When an exhaust jet calculation is made, the nozzle wall exit lip 
mesh point becomes a singularity and, therefore, is treated by a special 
procedure. This procedure consists of first computing an upstream solution 
at the exit mesh point using the flow tangency condition as the boundary 
condition and backward C differences in both the initial— value and solution 
planes. Next, a downstream solution is calculated using Eq. (13) as the 
boundary condition and the total conditions calculated from the upstream 
mesh point. The upstream solution is used when computing wall mesh points 
upstream of the exit mesh point, while the downstream solution is used 
when computing downstream wall mesh points. A third exit mesh point 
solution to be used for interior mesh point calculation is determined as 
follows. When the upstream solution is subsonic, the two solution Mach 
numbers are averaged such that the averaged Mach number is less than or 
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equal to one. This Mach number Is then used to calculate the exit mesh 
point solution to be used to compute the interior mesh points. When 
the upstream solution is supersonic t the upstream solution is used to 
calculate the interior mesh points. 

Stability 

The step size, At, is controlled by the well known Courant condition 
which can be expressed as 

At < l/[(V+a)(l/Ax^ + 1/Ay^)’"^] (14) 

where V is the velocity magnitude. Using Eqs. (5) and (10), Eq. (14) 
can be written as 

At < A/[(V+a)(l/Ac^ + (15) 

where the coefficient A was determined from actual calculations and varied 
between 1.0 and 1.6 depending on the geometry of the flow in question. 

Overall Program 

The nozzle inlet flow as well as the flow leaving the nozzle 
may be either subsonic or supersonic. The flow may contain stream-wise 
variations in stagnation temperature and stagnation pressure. The nozzle 
wall and centerbody geometries may be either one of three analytical con- 
tours or a completely general tabular contour. The program is capable of 
calculating the exhaust jet boundary for either subsonic or supersonic 
flow. The initial data may be either read in or calculated internally by 
the program. The internally computed data are calculated assuming one- 
dimensional, steady, isentropic flow with area change. The program output 
includes the coordinates, velocities, pressure, density, Mach number, 
temperature, mass flow, and axial thrust in both English and metric units. 
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Results and Discussion 


The results in the present study were obtained using a CDC 6600 com- 
puter. The computational times given are the central processor time not 
including compilation. In order to compare these results with those of 
other investigators, the following table of relative machine speeds is 
given. 

Computer Relative Machine Speed 


IBM 7094 

0.1 

IBM 360/50 

0.1 

IBM 360/65 

0.3 

IBM 360/75 

0,5 

Univac 1108 

0.5 

CDC 6600 

1.0 


These relative speeds were obtained from Refs. 13 and 14 and are only 
rough estimates. These values may vary considerably depending on the 
compiler, machine configuration, and numerical technique. The initial 
data in each case were the one-dimensional values computed internally by 
the program. When the relative change in axial velocity in the throat 
and downstream regions was less than a prescribed convergence tolerance, 
the flow was assumed to have reached steady state. The convergence 
tolerance was found to be a function of the mesh spacing, flow speed, and 
nozzle geometry. For the results presented here a convergence tolerance 
of 0.003 per cent for flows without exhaust jet calculations and 0.005 per 
cent for flows with exhaust jet calculations was employed. 

The present method was used to compute the steady state solution for 
flow in the 45® - 15® conical, converging-diverging nozzle shown in Fig. la. 
The Mach number contours and wall pressure ratio are shown in Fig, 2. 
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The experimental data are those of Cuff el, et al. The computed discharge 

coefficient is 0,983 as compared with the experimental value of 0.985. 

A 21 X 8 computational mesh was used, which required 301 time planes and 

a computational time of 35 seconds. In general, there is good agreement 

2 

with the experimental data. This case was also solved by Prozan , Migdal, 
et al. , Laval, and Serra. While the details of Prozan's computation were 
not reported by Cuffel, et al, , Saunders reported a time of 45 minutes on 
a CDC 3200 (23 x 11 mesh) for computing the flow in a nozzle with a large 
radius of curvature. Migdal, et al., reported a computational time of 
less than 5 minutes on an IBM 360/75. Laval reported a computational time 
on the order of 2 hours on an IBM 360/50 (61 x 21 mesh). Serra reported a 
computational time of 80 minutes on a Univac 1108 (3000 mesh points). In 
addition, this case was also solved by Prozan and Kooker^^ using a re- 
laxation scheme to solve the ir rotational equations of motion. They re 
ported a computational time of 5 to 10 minutes on an IBM 7094 (21 x 11 mesh). 

The present method was also used to compute the steady state flow in 
a 15® conical, converging nozzle. The nozzle geometry is shown in Fig. lb. 
The Mach number contours and wall pressure ratio for a nozzle pressure 
ratio of 2.0 are shown in Fig. 3. The experimental data are those of 
Thornock^^. The computed discharge coefficient is 0.957 as compared with 
the experimental value of 0.960. A 23 x 7 computational mesh was used, 
which required 249 time planes and a computational time of 29 seconds. 

In general, there is good agreement with the experimental data. This 
case was also solved by Wehofer and Moger and Brown and Ozcan. Wehofer 
and Moger* s solution for a pressure ratio of 2 required over 2 hours on 
an IBM 360/50 (47 x 11 mesh), while Brown and Ozcan's results required 17 
minutes on an IBM 360/65 (20 x 6 mesh). 
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Finally, the present method was used to calculate the flow in a 10® 
conical, plug nozzle. The nozzle geometry is shown in Fig. Ic. The 
Mach number contours and plug pressure ratio for a nozzle pressure 
ratio of 3.29 are shown in Fig. 4. The experimental data are those of 
Bresnahan and Johns^^. A 31 x 6 computational mesh was used, which required 
327 time planes and a computational time of 52 seconds. In general, there 
is good agreement with the experimental data. The author is unaware of 
any other time- dependent analysis of plug nozzles. 

Concluding Remarks 

A method of computing nozzle flows has been presented. A production 
type computer program capable of solving a wide variety of nozzle flows 
has been developed. The accuracy of this program was demonstrated by 
computing the flow in a 45® - 15® conical, converging- diverging nozzle, 
a 15® conical, converging nozzle, and a 10® conical, plug nozzle. The 
computational times were less than one minute on a CDC 6600 computer, 
which is considerably faster than any of the previous time-dependent 
techniques. 
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Figure 2, Mach number contours (above) and wall pressure ratio for 
AS'*- 15® conical noszle. 
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Figure 3, Mach nuni>er contours (above) and wall pressure ratio for 
15^ conical nozzle. 
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Figure 4. Mach number contours (above) and plug pressure ratio for 
10” conical, plug nozzle t 


